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We propose an experimentally relevant protocol for the controlled generation of matter-wave 

dark solitons in atomic Bose-Einstein condensates (BECs). In particular, using direct numerical 

simulations, we show that by switching-on a spatially inhomogeneous (step-like) change of the 

s-wave scattering length, it is possible to generate a controllable number of dark solitons in a quasi- 

C^ . one-dimensional BEC. A similar phenomenology is also found in the two-dimensional setting of 

(S| ■ "disk-shaped" BECs but, as the solitons are subject to the snaking instability, they decay into 

vortex structures. A detailed investigation of how the parameters involved affect the emergence and 

^ , evolution of solitons and vortices is provided. 
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INTRODUCTION 



c/:) ■ The experimental realization of Bose-Einstein condensates (BECs) of dilute atomic gases ^L] is one of the most 



fundamental developments in quantum and atomic physics over the last two decades. Many of the relevant theoretical 
and experimental studies [2|, |3| have been focused on the nonlinear excitations of BECs [J, |5|. Among them, the 
so-called matter-wave dark solitons (see the recent review [6|), namely localized density dips with a phase-jump across 
them, have been observed in BECs in a series of experiments |7l-[l7| . The interest on this type of localized nonlinear 
^ ■ structures may be attributed to a series of reasons: matter- wave dark solitons are one of the most fundamental 
Q^l nonlinear excitations of BECs with repulsive interatomic interactions, being the nonlinear analogs of the excited 
eigenstates of the quantum harmonic oscillator [18|, [l9|; furthermore, they arise spontaneously upon crossing the BEC 
2 \ phase-transition [20|, l21| , and their properties may be used as a diagnostic tool probing the rich physics of BECs [22| . 
Furthermore, matter- wave dark solitons may also be relevant to applications: for example, it has been proposed that 
the position of a dark soliton can be used to monitor the phase acquired in an atomic matter- wave interferometer in 
the nonlinear regime [23|, [2^ . 
Q ■ Matter- wave dark solitons have been generated in atomic BECs by means of various different methods. In particular, 

Q \ in the early experiments [3, Q (but also in some recent ones [Isl, llj), dark solitons were created by means of the 
phase-imprinting method: this technique relies on the illumination of part of the condensate by a short off-resonance 
laser beam, so as to imprint the necessary phase jump characterizing the dark soliton. Another technique, that was 
used in the experiments reported in Refs. [sl, [ll| (see also [25|) is the density- engineering method^ which is based on 
OO ' ^he creation of local density depletion on the BEC density by means of external potentials. A combination of the 
CO . phase-imprinting and density-engineering methods, namely the so-called quantum- state engineering technique is also 
^O ■ possible: this method involves manipulation of both the BEC density and phase, and has also been used in experiments 
^^ \ with two-component BECs [10|,[l3|- Furthermore, in Ref. [12|, the BEC flow against a broad penetrable barrier, which 
f-*^ ■ was swept through the condensate, led to the breakup of the BEC superfluidity and the concomitant generation of 
^D \ dark solitons for a particular regime of barrier speeds. Finally, in more recent experiments the the interference between 
^p ■ two (or more) condensate fragments confined in a trap was demonstrated to result in the formation of dark solitons 
|15l-[l7l, l25j , or vortices in higher-dimensional setups [26|, [g^ • 

In this work, we propose an alternative method to create dark solitons in quasi one-dimensional (ID) BECs, or 
vortices in higher- dimensional settings (particularly, in quasi two-dimensional (2D) BECs). The proposed technique 
relies on the manipulation of the s-wave scattering length, which controls the effective nonlinear coefficient in the 
J^ ■ Gross-Pitaevskii (GP) mean-field model: in particular, we consider a situation where the collisional dynamics across 
the condensate is inhomogeneous due to the presence of a spatially dependent external field in the vicinity of a 
Feshbach resonance. Note that such, so-called, collisionally inhomogeneous condensates 1281 , have attracted much 
attention, as they have been proved to provide a variety of interesting phenomena |29l-l41|. In the present study, 
we assume, more specifically, that by switching-on an external magnetic or optical field, the nonlinear coefficient 
becomes piecewise constant (in line with the considerations of [30i, i31l, SO, Ell ) , and can vary between two (smoothly 
connected) values. Then, depending on the magnitude of the jump in the value of the nonlinear coefficient, one or 
more dark solitons in the ID setting, or vortices in the 2D setting, can controllably be created. It is important to 
note that the considered form of the spatial change of the nonlinear coefficient (i.e., of the s-wave scattering length) 
may be straightforwardly experimentally realizable, through the use of magnetic field gradients of moderate size for 
atom chips (see also discussion in Ref. [41|). It should also be added that a similar technique has been proposed very 
recently in the three-dimensional context as a means of generating vortex rings [42|] . 



The paper is organized as follows. In Sec. |II]we present our model and discuss the dark soliton generation in the 
quasi-lD setting. In Sec. Illll we present results concerning the soliton generation, as well as their subsequent decay 
into vortices, in the quasi-2D setting and, finally, in SecjIVlwe present our conclusions. 

II. THE MODEL AND ITS CONSIDERATION IN THE ID SETTING 

A. Presentation of the model 

We consider a cigar-shaped (quasi- ID) condensate, confined in an anisotropic trap with longitudinal and transverse 
confining frequencies (denoted by uOx and o;^, respectively) such that cj^ <C o;^. In this case, taking also into regard 
the quasi-one-dimensional nature of the dynamics, it can be found [43|, Q ^^^^ ^^^ condensate evolution can be 
described by an effectively ID, GP-like equation which can be expressed in the following dimensionless form: 

idtu = —-dlu + V{x)u + \/l + g\u\'^u — /nu. (1) 

Here, u{x^ t) is the macroscopic wave function, normalized so that J_ \u\^dx = N (where N is the number of atoms), 
/i is the chemical potential, while the density \u\^ ^ the length, time and energy are measured, respectively, in units 
of (2a) ~^ (where a is the 5- wave scattering length), the transverse harmonic oscillator length a±_ = ^/h/muj±^ ^I^, 
and Jiw±. Note that the deviation from ID is accounted for by the generalized nonlinear ity in Eq. ([1]), which becomes 
the traditional cubic nonlinearity in the weakly-interacting limit of g\u\'^ <C 1 (in this limit, Eq. ([1]) is reduced to the 
usual ID GP equation - see, e.g., discussion in [4] and references therein). The potential V{x) in Eq. ^ is assumed 
to have the usual harmonic form, 

V{x) = \^^x\ (2) 

characterized by its strength Q. = uJxj^^- Finally, as concerns the nonlinear coefficient ^, it is considered to be of the 
following form: 



9{x',t) = {^4^ " ^^""^ ©1 }{M^ ^ ^^""^ (^7^)1 } ^ ^' ^^^ 



where the spatial and temporal scales s and r, the characteristic time to, and the constant £, are defined below, 
upon considering various asymptotic limits associated to the assumed form oi g{x\ t). In particular, in order to better 
understand the physical situation described by Eq. (|3|), let us first consider the asymptotic limit of t <C to: in this 
case, Eq. (J3]) describes a condensate with a constant nonlinear coefficient equal to ^ = 1. On the other hand, in the 
limit t ^ to, the nonlinearity coefficient asymptotes to the spatially dependent form ^o(^) = (1/2)^[1 — tanh(x/s)] + 1; 
this suggests, in turn, that for x -^ ±oo the nonlinearity coefficient takes the values g = 1 and g = l-\-e^ respectively. 
In other words, the condensate is initially (t <C to) characterized by a constant nonhnear coefficient ^ = 1, but 
for large times (t ^ to) the implementation of a temporal smooth ramp renders the nonlinear coefficient spatially 
inhomogeneous, experiencing a smooth jump, from ^ = l+£to^ = l, around the trap center {x = 0). According 
to the above, it is clear that e denotes the amplitude of the jump in the value of the nonlinear coefficient, s is the 
spatial length scale on which the transition between the values g = I and g = 1 -\- e takes place, to represents the 
characteristic time around which the spatial inhomogeneity of the nonlinear coefficient is "switched-on" , while the 
time scale r is the duration of this process. 

From a physical point of view, the above situation may be implemented experimentally as follows. First, the 
modification of the nonlinear coefficient g - i.e., of the scattering length a - from the value ^ = 1 to the value 
g = 1 -\- e can be achieved by using an external magnetic (or optical) field, which is realized around t = to- The 
desired inhomogeneity of the scattering length may be realized upon employing a bias homogeneous external field and 
imposing on top of it a steep localized spatial gradient. This leads to a constant scattering length in the left portion 
of the condensate (x < 0), followed by a localized change of a, ending up with a different value in the right portion 
of the condensate {x > 0). In our consideration, we choose the function tanh to mathematically characterize the 
switch-on of the external fields, as well as the transition between the different constant values of a. Also, it is natural 
to assume that we are relatively close to a Feshbach resonance, so as to easily manipulate the scattering length with 
the imposed external field. 
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FIG. 1: (Color online) Spatiotemporal contour plots plots, showing the evolution density |i^|^ of the condensate, as per Eq.([T]), 
with fixed chemical potential (/x = 2.5) and different values of the parameter e. When e = 0.8 dark solitons are not generated, 
while as the value of e is increased, the formation of 1, 2, 3, 4, 5 solitons is observed for s = 2, 3, 3.6, 4.2, 4.5, respectively. 



B. Dark soliton generation 



Let us now proceed with the presentation of our numerical results demonstrating the dark soliton generation in the 
ID setting considered above. The numerical simulations start with the determination of the stationary ground state 
of the system, obtained by means of Newton's method, with an initial condition corresponding to the value g = 1 
of the nonlinear coefficient. Then, the dynamical evolution of the system is monitored by a 4th-order Runge-Kutta 
method, with step At = 0.001. In the results presented below, the parameters are chosen as: ft = 0.02, to = 10, r = 1 
and s = 1. Notice that for the generation of dark solitons, the values of the chemical potential and of the jump of the 
nonlinearity coefficient (described by the parameters e and /i, respectively) are of crucial importance. 

In Fig. [H where some typical scenarios of the evolution are shown, we keep the value of the chemical potential fixed, 
namely /i = 2.5, and investigate how the value of e affects this process. We observe, at first, that dark solitons are 
not formed for sufficiently small values of e; a pertinent example, corresponding to £ = 0.8, is provided in the top- left 
panel of the figure. On the other hand, the increase of the value of e results in the formation of dark solitons. The 
top-middle panel of Fig. [T] shows that a dark soliton is emitted forwardly (i.e., toward the right portion of the BEG, 
with the smaller value of the scattering length) at about t = 14 when s = 2.0; similarly, the top-right panel of Fig. [1] 
shows the emergence of two solitons, for £ = 3, and so on (see a more detailed analysis below). 

Figure [1] also demonstrates the basic features of the condensate and soliton dynamics. In that regard, first we note 
that the switch-on (around t = to = 10) of the external potential modifying the spatial distribution of the nonlinear 
coefficient, results in the excitation of the quadrupole mode of the condensate [2|. In particular, after t = 10, the 
condensate sh ows a breathing behavior, performing the quadrupole oscillation with a characteristic frequency close to 
the value ~ y^5/2Q [4^. On the other hand, the dark solitons, when formed, perform oscillations due to the external 
potential. As seen in Fig. [H dark solitons with sufficiently large velocities are transmitted through the "interface" at 
the trap center (x = 0), while lower velocity solitons are chiefly reflected by the interface. Moreover, it can readily 
be observed that in cases where more than two solitons are formed, they may undergo elastic collisions: see, e.g., 
the third panel of Fig. [1] (for £ = 3), where the solitons collide at t ~ 22. Similar events arise in cases where more 
solitons are generated. It should also be mentioned that the soliton generation (and subsequent evolution) is followed 
by emission of radiation of weak linear (sound) waves, which are originally generated in the left portion of the BEG 
(with the larger value of the scattering length), but then travel in the whole system. 

At this point, it is necessary to mention that the emergence of solitons is a continuous process. In particular, 
instead of emerging only for a specific value of e, the density of a dark soliton becomes deeper as e increases and, 
meanwhile, the phase dislocation increases also, being associated with the presence of the dark soliton. Thus, if a 
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FIG. 2: (Color online) The number of solitons emitted before t = 25, according to the identification rule defined in section HIl 
as a step function of s for different values of the chemical potential /x. The five functions with various colors and line types 
represent five different cases chosen as example, i.e., /i = 2.0, 2.5, 3.0, 3.5 and 4.0. 

dark soliton is formed for a given value of £, the same soliton tends to appear earlier as e is increased. Due to the 
gradual character of this process, it is difficult to quantitatively identify the emergence of the soliton(s). Nevertheless, 
our rule of identification is set in the way that, if a dark soliton is emitted forwardly {x > 0) before t = 25, then it is 
counted. This is how we determine the number of solitons as in Fig. [H 

As mentioned before, besides £, the chemical potential /i is another factor which is important for the soliton 
formation. Figure [2] demonstrates the joint influence of these two parameters. The increase of each parameter has a 
positive effect on the system's capability of emitting solitons. In particular, the larger the value of /i is, the smaller 
the value of e at which dark solitons may be generated. For instance, applying the identification rule stated before, 
the second soliton emerges when e = 3.08 for /i = 2.0, or when e = 2.41 for /i = 3.0, and so on (see Fig. [2j). 



III. THE 2D SETTING 



We have also studied the 2D setting, corresponding to the case of the so-called disk-shaped BECs. In this case, the 
evolution of the BEC will be described by the pertinent GP equation with the cubic nonlinearity [2|, Q : 



idfU 



--Au + V{r)u + g\u\'^u — fiu^ 

2^2 



where u{x^y^t) represents the wave function, V{r) = (l/2)Vt r is the harmonic potential with r 



(4) 

x^ + ^^, and 
the Laplacian A = 9^ + 9^ [48|. The nonlinear coefficient g is assumed to take the same form, as in the ID setting, 
indicating that the spatial inhomogeneity is employed along the x-axis. 

As before, we start the simulations with the steady state of the system, corresponding to ^ = 1. The parameters s, 
to and r remain the same as in the ID setting. During the evolution, dark solitons are emitted when the amplitude e 
and chemical potential fi are sufficiently large. Moreover, another parameter is found to play an important role in this 
process, namely the normalized strength Vt of the trapping potential V{r). In the numerical results of this section, ji is 
fixed to 3.0, as its influence on the dark soliton has been discussed in detail in the ID case, and we focus on cases with 
various pairs of Q. and e. Figure [3] describes a typical evolution scenario, by showing snapshots of the BEC density 
for Vt = 0.08 and e = 2.5. Starting with a stationary ground state (similarly to the ID case), the condensate becomes 
anisotropic due to the change of the value of the nonlinear coefficient; a rectilinear dark soliton is generated near the 
line X = 0, and then it moves to the forward direction (x > 0), as shown in the top panels of the figure. Subsequently, 
this rectilinear dark soliton does not oscillate as it does in the ID setting; instead, it undergoes a bending and, after 
moving forward for a short time, it finally breaks up into a sequence of vortex pairs, as a result of the onset of the 
snaking instability (see, e.g., [6] and references therein, as well as [10| for relevant experimental observations) - see 
the middle row of Fig. [3l The generated vortices are moving due to their interactions and the action of the external 
potential, and some of them disappear at later times. It is also worth observing (see the bottom row of Fig. [3]) the 
manifestation of an interesting transient phenomenon, namely the formation of a closed-loop structures, reminiscent 
of ring dark solitons in BECs [45] which, however, eventually disappear. After a certain time, no specific pattern is 
captured in the density plot, as shown in the last panel of the figure. 

In order to monitor the overall evolution of the vortices in the system, we compute the vorticity based on the 




FIG. 3: (Color online) Contour plots showing snapshots of the evolution of the 2D condensate density for Q 
(the chemical potential is /i = 3). 



: 0.08 and s = 2.5 



corresponding fluid velocity (see, e.g., Ref. [46[), namely. 



ii*Vii - iiVii* 



Ve = -Z- 



The fluid vorticity is then defined as 



u[ 



V X Ve 



(5) 



(6) 



According to the above, the evolution of the vortex structures depicted in Fig. [3] is presented in the top right panel of 
Fig.m The accompanying top left panel of the same figure is a spatiotemporal contour plot showing the evolution of the 
condensate density along the x-direction (i.e., on the plane ?/ = 0), while the middle and bottom panels correspond 
to different (larger) values of the parameters Q and e. This figure provides a view complementary to the one of 
Fig. [3] concerning the emergence of the dark soliton and vortices, and their subsequent evolution. It is interesting to 
highlight that while, as indicated above, no specific density pattern seems to be selected by the asymptotic dynamics, 
nevertheless, it is clear that the topological charge (ultimately) induced by the nonlinearity step clearly seems to 
persist throughout the timescales that monitored herein. 

For a given Q, more dark solitons are emitted when e increases. As expected, and similarly to the above case, the 
solitons decay again to vortices due to the onset of the snaking instability: in fact, this happens shortly after their 
emergence (and their subsequent forward motion for short time). Here, it should be mentioned, however, that the 
evolution of the vortices differs in each different case. To be more specific, we provide two more examples, one with 
three solitons and the other with two solitons, shown in the middle and bottom panels of Fig. |4] (parameter values are 
^ = 0.1 and 0.12, e = 5 and 7.5, respectively). A series of snapshots of the case corresponding to the middle panels 
of Fig. m is provided in Fig. [5l which demonstrates the emission of the two dark solitons, together with the evolution 
of vortices. 



IV. CONCLUSIONS 



In this work, we have proposed and analyzed an experimentally relevant technique for the creation of dark solitons 
and vortex structures in Bose-Einstein condensates. Our method, can be briefly described as follows. After creating a 
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FIG. 4: (Color online) Contour plots showing the evolution of the density on the y = plane (left panels) and the spatiotemporal 
contour plots of the vorticity (right panels) - see text. From top to bottom, parameter values are Q = 0.08, 0.1 and 0.12, 
e — 2.5, 5.0 and 7.5, respectively (the chemical potential is /x = 3). 



condensate, characterized by a specific scattering length, an external magnetic or optical field (with a value relatively 
close to a Feshbach resonance) is switched-on. If this external field is spatially inhomogeneous, i.e., it has a steep 
localized spatial gradient on top of a constant value, then the scattering length becomes spatially inhomogeneous too, 
taking different values to the left and to the right of the trap center. We have used a mean- field model in (1 + 1)- 
and (2 + l)-dimensions, describing cigar-shaped and disk-shaped condensates, respectively, in order to investigate the 
effect of such a collisional inhomogeneity. 

Using direct numerical simulations, we have shown that if the change of the value of the scattering length exceeds 
a certain threshold, matter-wave dark solitons are spontaneously generated in the system by means of a continuous 
process. We have found how this threshold for soliton generation depends on other parameter values, e.g., in the 
quasi- ID setting it becomes smaller when the chemical potential (i.e., the number of atoms) is increased. We have 
also shown that rectilinear dark solitons can also be spontaneously formed in a disk-shaped condensate. Nevertheless, 
due to the higher dimensionality of this quasi-2D condensate, the dark solitons are subject to the snaking instability. 
Thus, shortly after their formation, they decay into vortex pairs. In certain cases, we have also observed the formation 
of interesting transient structures, in the form of closed-loop patterns, reminiscent of ring dark solitons. 

At this point, it is relevant to mention that, very recently, a detailed experimental control of the inter- atomic 
interactions in a Bose-Einstein condensate via optical Feshbach resonances was reported [43. We believe that, as 
further experiments in the same direction are expected to appear, our results are likely to be useful as a systematic 
means of producing fundamental nonlinear excitations in such collisionally inhomogeneous atomic systems. 





FIG. 5: (Color online) Evolution of the 2D condensate density with Q = 0.1 and e = 5.0 
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